A Dynamic Interstellar Medium: 
Recent Numerical Simulations 

Anvar Shukurov 

Department of Mathematics, University of Newcastle, 
Newcastle upon Tyne NE1 7RU, England 

Abstract. Recent numerical simulations of the interstellar medium driven by energy input from 
supernovae and stellar winds indicate that H I clouds can be formed by compression in shock waves 
and colliding turbulent streams without any help from thermal and gravitational instabilities. The 
filling factor of the hot phase in these models does not exceed 20—40% at the midplane. The hot gas 
is involved in a systematic vertical outflow at \z\ < 1-3 kpc, similar to that expected for galactic 
fountains, whereas the warm component may remain in hydrostatic equilibrium. The turbulent ve- 
locity is larger in the warmer phases, being 3, 10 and 40kms _1 in the cool, warm and hot phases, 
respectively, according to one of the simulations. The models exhibit global variability in the total 
kinetic and thermal energy and star formation rate at periods of (0.4—4) X 10 8 yr. Current models 
are still unable to reproduce dynamo action in the interstellar gas; we briefly discuss implications of 
the dynamo theory for turbulent interstellar magnetic fields. 



1 Introduction 

Recent numerical simulations of the interstellar medium (ISM) heated by supernovae 
and stellar winds have highlighted the importance of numerical experiment in the 
studies of the multiphase ISM. The simulations play a complimentary role to obser- 
vational and semi-analytical studies and are especially important in clarifying such 
long-standing problems as the filling factors of the various phases of the ISM, the 
nature and lifetimes of H i clouds, dynamo action, etc. Simulations capturing relevant 
physical effects at appropriate resolution have become feasible only recently, and even 
first results reviewed here provide insight into many controversial issues. 

The range of gas temperatures and densities in the ISM is tantalizingly wide, 
from T ~ 10 6 K, n ~ 10~ 3 cm -3 in the hot phase to 10 K and more than 10 3 cm -3 in 
molecular clouds. The gas is involved in random motions over a wide range of scales 
(e.g., Spangler, 1999). It is still impossible to cover the whole range of parameter 
variation (i.e., all the phases) in a single model of the ISM. In this paper I shall 
review models of the diffuse phases of the ISM with densities below 10-100 cm -3 . 
Simulations of MHD turbulence in dense clouds have been discussed by Padoan et al. 
(1997, 1998), Stone et al. (1998), Ostriker et al. (1999) and Heitsch et al. (1999); an 
extensive review and references can be found in Vazquez-Semadeni et al. (2000). 
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2 Models of the multiphase ISM driven by supernovae 

In this paper we discuss the results of three models of the ISM briefly summarized 
in Table 1. The models include similar physical effects and their results agree in 
many respects. All are non-linear, hydrodynamic or magnetohydrodynamic models 
that include external gravity (typical of the Solar neighbourhood), suitably idealized 
heating sources and radiative cooling (assuming partially ionized, optically thin gas). 
All models employ artificial viscosities to avoid excessively large gradients in the 
simulated quantities. Neither of the models includes explicitly ionization balance and 
cosmic rays. 



Tab. 1. Recent comprehensive numerical models of the ISM driven by stellar winds and supernovae. 

RBN (a) 



Dimension 
Box size, kpc 
Resolution, pc 
Heating 

Temperature range, K 
Magnetic fields 
Self-gravity 
Rotation 
Cosmic Rays 
Ionization Balance 



2D: xy and xz 
2 x 2 and 2 x 15 
10 and 10 x (10-1000) 
SN II, winds 
10 2 -10 8 



VSPPW 
2D: xy 
1 x 1 
2 

Diffuse, winds 
10 2 -4 x 10 4 
imposed 
yes 
yes 



KBSTN (c) 
3D: xyz 
0.5 x 0.5 x 2 



SN II, SN I 

10M0 8 
self-excited 

yes 



References 

( a ) Rosen et al. (1993, 1996), Rosen & Bregman (1995); 

( b ) Vazquez- Semadeni et al. (1995), Passot ct al. (1995), Ballestcros-Paredes et al. 

(1999); 
MKorpi et al. (1999a,b). 



The model of Rosen et al. (1993, 1996) and Rosen & Bregman (1995) (labelled 
RBN) uses a finite- difference ZEUS code in two dimensions and neglects magnetic 
fields and rotation. These authors have results for a region in the Galactic plane (xy) 
and in a vertical plane (xz); the mesh size in the z-coordinate is variable, growing 
with z from 10 pc to 1 kpc. A unique feature of this model is that the stellar Pop- 
ulation I is included explicitly via equations for the stellar fluid coupled to the gas. 
The heating sources are Type II SNe, which inject energy instantaneously, and stellar 
winds modelled as a continuous heat source; both are correlated with the stellar den- 
sity. To alleviate numerical problems, the simulations have been split into alternating 
stages, the interaction phase of the gas and stars when any motions arc neglected, 
and a stage of hydrodynamic evolution when any coupling between stars and gas is 
neglected. SN explosions produce very hot gas, so the model contains three phases, 
i.e., the hot, warm and cool gas. 

Vazquez-Semadeni et al. (1995, 1996, 1997, 1998), Passot et al. (1995), Scalo et 
al. (1998) and Ballesteros-Paredes et al. (1999), labelled VSPP in Table 1 (sec also 
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earlier results in Passot et al., 1988, Lcorat ct al., 1990, and Vazquez-Semadeni & 
Gazol, 1995), employ a pseudo-spectral numerical scheme and restrict themselves to 
two dimensions in the Galactic plane. Their best runs have a spatial resolution of 2 pc 
or less, but only features at scales above ~ 10 pc are reliable (Passot et al., 1995). 
To avoid numerical difficulties, this model includes diffuse heating (which is thought 
to model the UV heating of the ISM) and stellar heating, but each stellar energy 
release event is spread over a period of 6 x 10 6 yr. Furthermore, both the heating and 
cooling rates are reduced by a factor 7-10 below the realistic values. As a result, only 
the warm and cool gas phases are generated, with the maximum gas temperature of 
4 x 10 4 K. A peculiar feature of this model is the inclusion of self-gravity; however, its 
role is only marginal at the densities reached in the simulations. The model includes 
magnetic field, but any dynamo action is precluded in two dimensions. 

The model of Korpi et al. (1999a, b) (KBSTN) is fully three-dimensional and in- 
cludes the effects of rotation and magnetic field at a resolution of 8 pc. This is the 
only of the three models admitting realistic behaviour of vorticity and magnetic field 
including vortex stretching and dynamo action. The gas is heated by SNe modelled 
as thermal energy release over one time step which can be as short as 10-100 yr. As 
a result, the gas is heated to T ~ 10 8 K at the explosion site and hot, warm and 
cool gas phases can be identified in the simulations. Both Type I and II SNe are 
included, which differ in both the occurrence rate and the vertical distribution. The 
SN explosions occur at randomly chosen sites, but Type I SNe can only explode in 
those regions where the gas density exceeds the average at that height. With this 
prescription, about 70% of the Type I SNe are clustered resembling OB associations. 
The simulations have only been run over time span of the order of 10 8 yr, so it is yet 
unclear whether or not any magnetic fields can be supported by dynamo action in this 
model. A related model is presented by Gudiksen (1999) where an intricate algorithm 
for choosing SN explosion sites involving modelling of the stellar initial mass function 
has been developed. 

Despite different numerical approaches employed and certain differences in the 
physical content, these models yield many coherent results which appear to be model- 
independent and which are our subject in what follows. 

3 The multi-phase structure and cool gas filaments 

In all the models, the cooling function is truncated at T = 100-300 K to avoid thermal 
instability at T < 10 5 K. Self-gravity is either neglected or unimportant. Nevertheless, 
cool, dense gas clouds with T < 10 3 K, n = 1-100 cm -3 are a typical feature of 
all the simulations. The clouds are elongated and are better described as filaments 
or sheets (which can be distinguished from filaments only in 3D models); rounder 
structures occur at the intersections of the filaments. The filamentary Hi clouds occur 
at positions of converging turbulent flow, V ■ v < and are not gravitationally bound 
(Vazquez-Semadeni et al., 1995; Ballesteros-Paredes et al., 1999). This indicates that 
they are formed by compression (ram pressure) in the turbulent velocity field. The 
lifetime of the filaments is shorter than 10 8 yr (Rosen & Bregman, 1995), which is close 
to the kinematic time scale. The filaments are mainly destroyed by mutual collisions. 
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In models with SNe, a more obvious mechanism of the cloud formation is compression 
in expanding shocks (Korpi et al., 1999). The length of the filaments is 50-100pc and 
they are often extended vertically because of the vertical outflow of the hot gas, 
resembling Hi 'worms' observed in spiral galaxies (Rosen & Bregman, 1995; Korpi et 
al., 1999). The three-dimensional model of Korpi et al. (1999) has both filaments and 
sheets of cool gas. In the model of Rosen & Bregman (1995), where SNe are confined 
to a rather thin layer, the cool filaments are not uncommon at z = 1-3 kpc, but in 
the model of Korpi et al. (1999) they are confined to \z\ < 100 pc because Type II 
SNe that can occur at large heights prevent cool gas from rising higher. 

As could be expected, the hot phase (T > 10 5 K) can only be generated by the 
SNe. The hot gas fills bubbles surrounded by dense shells; it breaks through the warm 
layer and streams into the halo (or beyond the upper boundary of the computational 
domain) at systematic velocities of 100-200 kms -1 (Korpi et al., 1999); Rosen & 
Bregman (1995) obtain 400 kms" 1 at z = 2-3 kpc. 

The volume filling factors of the hot, warm and cold phases are 40-30%, 40-60% 
and 20-10% at z = 1-3 kpc in the model of Rosen & Bregman (1995). The filling 
factor of the hot gas increases at larger heights at the expense of the other phases. 
The filling factor of the hot phase in the model of Korpi et al. (1999) grows from 
20-30% at z = to 80-100% at \z\ = lkpc; the remaining part of the volume is 
mostly occupied by the warm gas. Unlike the hot gas, the warm component is, on 
average, in hydrostatic equilibrium with a scale height of 200 pc; this agrees with 
a detailed analysis of hydrostatic equilibrium in the Solar vicinity by Fletcher & 
Shukurov (1999). As mentioned above, the cold gas is confined to \z\ < 100 pc. Rosen 
& Bregman (1995) perform detailed fitting of the vertical gas distribution in different 
phases. The scale heights of the cold, warm and hot gas in one of their models (Run E) 
are 225, 550 and 2000 pc, respectively. Note that the gas called 'cold' in the model, 
can be identified with the cold and warm neutral medium, whereas the 'warm' gas 
corresponds to the warm ionized medium. 

The phases of the modelled ISM are in a rough thermal pressure equilibrium. For 
example, density and temperature span four orders of magnitude in the model of 
Korpi et al. (1999), but thermal pressure varies in space only by an order of magni- 
tude. Similarly, the density contrast is a factor of 50 whereas the pressure contrast is 
only a factor of 5 in the simulations of Vazquez-Semadeni et al. (1995). The approx- 
imate pressure equilibrium is not trivial to explain since the sound crossing time is 
not shorter than the dynamic time scale because both the turbulent and systematic 
motions are transonic or supersonic. Vazquez-Semadeni et al. (1995) discuss how the 
pressure balance can result from the instantaneous establishment of thermal balance 
in the gas. With their diffuse heating, both the cooling time and the heating time 
are 10-100 times shorter than the sound crossing time and the dynamic time scale, 
10pc/10kms _1 ~ 10 7 yr. Therefore, the heating rate T must be always close to the 
cooling rate pA (with p the gas density) ; for A oc T m (with m a temperature-dependent 
index), the equilibrium temperature follows as T eq oc (T/p) 1 /™ 1 . The resulting pressure 
(assuming ideal gas) is P cq cx T 1 / m p y , where 7 = 1 — 1/m = 0.33-0.66 for m = 1.5-2.9. 
Since 7 < 1, denser regions are cooler, and this results in weaker spatial variations of 
pressure (Vazquez-Semadeni, 1998; Vazquez-Semadeni et al., 2000). This explanation 
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relies on a non-localized nature of the heating (it also applies if the heating rate is 
weakly dependent on density — Passot et al., 1995), but it may be of more general 
importance. 

4 Global variability 

Quasi-periodic variations in the total thermal and kinetic energy (and also magnetic 
and gravitational energy where appropriate) accompanied by similar variations in star 
formation intensity are typical of the models. The period of the variations changes 
from one model to another in the range (0.4-4) x 10 8 yr. Rosen & Bregman, 1995, do 
not discuss the global variability but mention oscillations, at a period 2.5 x 10 8 yr, of 
the extent of the region occupied by filamentary gas structures. This is compatible 
with the vertical sound crossing time and free-fall time across 1 kpc, so these may 
be a global acoustic mode, although Vazquez-Semadeni et al. (1995) interpret the 
variations as an analogue of spiral density waves in their model. 

5 Turbulence in the multiphase ISM 

The models include a random element in the heating source, so it is not surprising that 
the gas motion is random, resembling developed turbulence. Vazquez-Semadeni et al. 
(1995, 1997) show that compressible motions with a kinetic energy spectrum Ek oc fc~ 2 
characteristic of shock waves contribute significantly to the motions. Explosive heating 
by SNe driving shock waves would further increase this contribution. These authors 
argue that the velocity dispersion-size relation for interstellar clouds, Aw oc R 1 / 2 may 
result not from any kind of virial equilibrium, but rather directly from the turbulent 
spectrum: v\ ~ kE^ oc where is velocity at a scale R oc k^ 1 . 

The difference between cloud formation mechanisms in expanding SN shells and in 
colliding turbulent streams becomes significant when the vorticity of the motions, w, 
is considered. Cool clouds formed by obliquely colliding streams of the turbulent flow 
must have enhanced vorticity: Vazquez-Semadeni et al. (1995) obtain an enhancement 
by a factor of 10 for denser clouds. Such an enhancement is not observed in a 3D 
model of Korpi et al. (1999a, b) where gas filaments and sheets can be also formed 
by expanding SN remnants and where essentially three-dimensional mechanisms are 
efficient in amplifying vorticity in all the phases. Vorticity is amplified at curved shock 
fronts distorted by ambient density inhomogencitics, by vortex stretching and by the 
baroclinic effect, and 60-90% of the turbulent energy is in vortical motions (Korpi et 
al., 1999a). Vortex stretching is an essentially three-dimensional effect; furthermore, 
the tendency for conservation of ui/ p in a two-dimensional flow also contributes to an 
artificial enhancement of vorticity in denser regions of a 2D model (Vazquez-Semadeni 
et al., 1995). 

Korpi et al. (1999) computed the autocorrelation function of the vertical veloc- 
ity for the cold, warm and hot phases separately. The total (three-dimensional) r.m.s. 
turbulent velocity of the cool gas is as low as 3kms _1 , whereas it is close to lOkms -1 
for the warm gas (approximately the speed of sound in the warm phase) . The corre- 
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lation radius in the warm gas is about 30 pc (this is half the eddy size) and it only 
weakly changes with z; this is interpreted as an indication of approximate vertical 
hydrostatic equilibrium of the warm gas with a scale height of about 200 pc. The 
r.m.s. random velocity in the hot gas, 40kms -1 at \z\ < 1 kpc, is significantly larger 
than in the other phases, but still remains significantly smaller than the sound speed 
in the hot phase, lOOkms -1 . The vertical streaming of the hot gas discussed above 
can feed turbulence higher in the halo, so the turbulent velocity of 60 km s _1 observed 
in the halo by Kalberla et al. (1998; see also Pietz et al., 1998) seems to be perfectly 
compatible with the model. (Turbulent velocities of order lOOkms -1 in the halo have 
been invoked in galactic dynamo models involving the halo — Sokoloff & Shukurov, 
1990; Brandenburg et al., 1992, 1993.) The typical horizontal radius of the hot re- 
gions at the midplanc is about 20 pc; this can be identified with a chimney radius. 
The correlation radius for the hot phase grows with z together with the volume filling 
factor. 

6 Magnetic fields and dynamo action 

The only model admitting hydromagnetic dynamo action among the three models 
discussed is that of Korpi et al. (1999a, b). These simulations are initialized with a 
weak (0.1 /iG) regular magnetic field to provide a seed field for the dynamo. Both the 
large-scale and the random magnetic fields show a tendency to grow in the model, 
but runs spanning at least 10 9 yr are required to demonstrate any dynamo action 
convincingly. 

Magnetic fields in the 2D simulations of Passot et al. (1995) do not decay only 
because the model has an imposed uniform magnetic field, maintained throughout the 
simulations, which is tangled by small-scale motions. Magnetic fields in the cool fila- 
ments formed by converging flows often have a reversal along the filament axis. Mag- 
netic energy in these simulations is close to the kinetic energy in solenoidal motions at 
all scales. Magnetic fields, like vorticity, are especially sensitive to the dimensionality 
of the model, and two-dimensional results should be considered with caution. 

According to dynamo theory, the growth time of the random magnetic field due to 
the small-scale dynamoH action is expected to be as short as the eddy turnover time, 
~ 10 7 yr in the warm phase. A plausible statistically steady state resulting from the 
saturation of the dynamo action is an ensemble of magnetic flux ropes (e.g., Zeldovich 
et al., 1983, Sect. 8. IV). The length of the rope is of the order of the turbulent scale, 

— 1/2 

I ~ 50-100 pc, and the thickness is of the order of lR m , where R m is the magnetic 
Reynolds number. Dynamo action can occur provided R m > i? mjC r, where the critical 
magnetic Reynolds number is estimated as i? m ,cr = 60-100 in simplified models of 
homogeneous, incompressible turbulence. 

Subramanian (1999) suggested that a steady state, reached via the back-action of 
the magnetic field on the flow, can be established by the reduction of the effective 
magnetic Reynolds number R miC ff down to the value critical for the dynamo action. 

1 In this context, 'small' scales are understood as those smaller than the correlation scale of the 
turbulence, 50— 100 pc in the warm gas. 
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Then i? mj0 ff « i? m ,cr in the steady state. Therefore, the thickness of the ropes in the 

— 1/2 

steady state can be estimated as lR m ,cr ■ Using a model nonlinearity in the induction 
equation with incompressible velocity field, Subramanian (1999) showed that the mag- 
netic field strength within the ropes bo saturates at the equipartition level with kinetic 
energy density, 6q/87t ~ ^pv 2 , where p is the gas density and v is the r.m.s. turbulent 
velocity. The average magnetic energy density is estimated as (b 2 /8ir) ~ CR^ cr ^pv 2 , 
where C is a numerical coefficient, presumably of order unity. The volume filling fac- 
tor of the ropes then results as / ~ R^ CI — 0.01; correspondingly, the mean magnetic 
energy generated by the small-scale dynamo in the steady state is about 1% of the 
turbulent kinetic energy density. 

Using parameters typical of the warm phase of the ISM, this theory predicts that 
the small-scale dynamo would produce magnetic flux ropes of the length (or the 

— 1/2 

curvature radius) of about I ~ 50-100 pc and thickness lR m ,c r ~ 5-10 pc. The field 
strength within the ropes, if at equipartition with the turbulent energy, has to be 
of order 1.5 /zG in the warm phase (n = 0.1cm -3 , v = lOkms -1 ) and 0.5/iG in 
the hot gas (n = 10~ 3 cm~ 3 , v = 40kms~ 1 ). Note that other models of the small- 
scale dynamo admit solutions with magnetic field strength within the ropes being 
significantly above the equipartition level (Belyanin et al., 1993). 

The small-scale dynamo is not the only mechanism producing turbulent magnetic 
fields (e.g., Beck et al., 1996, Sect. 4.1 and references therein). Any mean-field dy- 
namo action producing magnetic fields at scales exceeding the turbulent scale also 
generates small-scale, turbulent magnetic fields. Similarly to the mean magnetic field, 
this component of the turbulent field presumably has a filling factor close to unity in 
the warm gas and its strength is expected to be close to equipartition with the turbu- 
lent energy at all scales. (As argued below, this component of the turbulent magnetic 
field may be confined to the warm gas, so magnetic field in the hot phase may have 
a better pronounced ropy structure.) 

The overall structure of the interstellar turbulent magnetic field in the warm gas 
can be envisaged as a quasi-uniform fluctuating background with one percent of the 
volume occupied by flux ropes (filaments) of a length 50-100 pc with well-ordered mag- 
netic field. The ropes can plausibly produce elongated, localized filamentary structures 
of polarized emission and Faraday rotation recently observed in the Milky Way by 
Wieringa et al. (1993) and Gray et al. (1998). The basic distribution described would 
be further complicated by compressibility, shock waves, MHD instabilities (such as 
Parker instability), etc. 

The site of the mean-field dynamo action is plausibly the warm phase rather than 
the other phases of the ISM. The warm gas has a large filling factor (so it can occupy a 
simply connected global region), it is, on average, in a state of hydrostatic equilibrium, 
so it is an ideal site for both the small-scale and mean-field dynamo action. Molecular 
clouds and dense H i clouds have too small a filling factor to be of global importance. 
Fletcher & Shukurov (1999) argue that, globally, molecular clouds can be only weakly 
coupled to the magnetic field in the diffuse gas. (However, the small-scale dynamo 
can work within the clouds.) The time scale of the small-scale dynamo in the hot 
phase is > l/v ~ 10 6 yr for I = 40 pc and v = 40kms~ 1 . This can be shorter 
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than the advection time due to the vertical streaming, h/V z ~ 10 yr with h = 1 kpc 
and V z = 100 km s" 1 . Therefore, the small-scale dynamo action should be possible in 
the hot gas even at small heights. However, the growth time of the mean magnetic 
field must be significantly longer than l/v reaching a few hundred Myr (e.g., Beck et 
al., 1996, Sect. 4.4). Thus, the hot gas can hardly contribute significantly to the the 
mean-field dynamo action in the disc and can drive the dynamo only in the halo. The 
fountain flow rather pumps the mean magnetic field out from the disc (Brandenburg 
et al., 1995); this effect can contribute to the saturation of the mean- field dynamo in 
the disc. 
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